home *** CD-ROM | disk | FTP | other *** search
/ Graphics Plus / Graphics Plus.iso / general / fractal / kaos.lha / fixptlib / fmapi_user.c < prev    next >
Encoding:
C/C++ Source or Header  |  1990-01-27  |  2.7 KB  |  108 lines
  1. /*
  2. ### Compute an inverse of f ^ r with an offset "xoffset" by secant method ###
  3. (Solve F(x)=f^r(x)-x-xoffset = x0)
  4. ------------------------------------------------------
  5. Note:    x[0...n-1] convention 
  6.     All computation in Euclidean coords
  7.     
  8. Input:    x[n] = current guess, n= phase space dimension, r= period
  9.     x0[n] = current point 
  10.     tolx= tolerated error in x (sum), tolf= tolerated error in F(x) (sum)
  11.     ntrial = # of maximum secant step
  12. Output:    x[n]=new guessed solution, *ndid= acutal # of secant step done,
  13.     *err= actual error in x (sum)
  14. */
  15.  
  16. #define DELTAX 1.e-6
  17. #define DELFRAC 0.2
  18. #define FREERETURN {*ndid = k;free_dmatrix(alpha,0,n-1,0,n-1);free_dvector(xt,0,n-1);free_dvector(beta,0,n-1);\
  19.     free_ivector(indx,0,n-1);return;}
  20.  
  21. void fmapi_user(x,x2,x0,ntrial,n,r,tolx,tolf,err,ndid)
  22. int ntrial,n,r,*ndid;
  23. double *x,*x2,*x0,tolx,tolf,*err;
  24. {
  25.     int i,j,k,*indx,*ivector();
  26.     double tolx10,fabs(),errf,d,*xt,*beta,**alpha,*dvector(),**dmatrix();
  27.     void usrfuni(),ludcmp(),lubksb(),free_dmatrix(),free_ivector(),free_dvector();
  28.     extern int stop,cur_color;
  29.     extern int model,mapping_on,fderiv_on;
  30.     extern double cutoff,*t_vf;
  31.  
  32.     tolx10 = tolx /10;
  33.     indx = ivector(0,n-1);
  34.     xt = dvector(0,n-1);
  35.     beta = dvector(0,n-1);
  36.     alpha = dmatrix(0,n-1,0,n-1);
  37.  
  38.     for(k=0;k<ntrial;k++){
  39.         for(i=0;i<n;i++){
  40.             if(x[i] > cutoff || x[i] < -cutoff) {
  41.  
  42.                 system_mess_proc(1,"inverse appears to blow up! Stop!");
  43.                 stop = 1;
  44.                 FREERETURN;
  45.             }
  46.         }
  47.         /* compute Jacobian by finite difference only */
  48.         (void) usrfuni(x,x2,x0,alpha,beta,r,n);
  49.         /*
  50.         printf("NTRIAL=%d\n",k);
  51.         printf("usrfun x: ");
  52.         for(i=0;i<n;i++) printf("%g ",x[i]);    
  53.         printf("\n");
  54.         printf("usrfun f: ");
  55.         for(i=0;i<n;i++) printf("%g ",beta[i]);    
  56.         printf("\n");
  57.         printf("usrfun alpha:");
  58.         for(j=0;j<n;j++){
  59.             for(i=0;i<n;i++) printf("%g ",alpha[j][i]);    
  60.             printf("\n");
  61.         }
  62.         */
  63.  
  64.         errf = 0.0;
  65.         for(i=0;i<n;i++) errf += fabs(beta[i]);
  66.         if(errf <= tolf) FREERETURN
  67.  
  68.         ludcmp(alpha,n,indx,&d);
  69.         /* should be placed here not to proceed further
  70.         in case of singular matrix (singular matrix often
  71.         arise when all the machine accuracy was lost in
  72.         computing alpha. This happens in case of good
  73.         convergence,too) */
  74.         if(stop){
  75.             /* do not reset stop=0 since it is passed down to
  76.             another routine */
  77.             FREERETURN
  78.         }
  79.         /*
  80.         for(j=0;j<n;j++){
  81.             for(i=0;i<n;i++) printf("%g ",alpha[j][i]);    
  82.             printf("\n");
  83.         }
  84.         */
  85.         lubksb(alpha,n,indx,beta);
  86.  
  87.         /*
  88.         printf("new beta:");
  89.         for(i=0;i<n;i++) printf("%g ",beta[i]);    
  90.         printf("\n");
  91.         */
  92.         *err = 0.0;
  93.         for(i=0;i<n;i++){
  94.             *err += fabs(beta[i]);
  95.             x2[i] = x[i];
  96.             x[i] -= beta[i];
  97.         }
  98.         if(*err <= tolx) FREERETURN
  99.         for(i=0;i<n;i++) if(fabs(beta[i]) < tolx10) x[i] -= tolx10;
  100.  
  101.         /* draw intermediate steps */
  102.         /*
  103.         (void) draw_record_pwf(x,cur_color,1,3,1,0);
  104.         */
  105.     }
  106.     FREERETURN
  107. }
  108.